本文为学习Möller-Trumbore算法的笔记，这里不对前面重心坐标求交进行记录，较为简单可自行推导。

根据中心坐标系可以很轻松的得出

$P = (1 - u - v)A + uB + vC$

这里P为射线与三角形相交点。u, v为以三角形两边为坐标系轴的P的投影长度(也可以理解为是P与该边两点连线三角形的面积与整体之比)。

$P=A - uA - vA + uB + vC = A + u(B - A) + v(C - A)$

而$P=O+tD$，因此代入可得

$\begin{array}{l}$
O+tD & = & A + u(B - A) + v(C - A)\\
O-A & = & -tD+u(B-A)+v(C-A)
\end{array}

这里转换为线性方程组的形式

$\begin{bmatrix}$
-D & (B-A) & (C-A)
\end{bmatrix}
\begin{bmatrix}
t\\u\\v
\end{bmatrix}
=O-A

那么问题就转化为如何求解出这里的$t$ $u$ $v$三个分量

运用克莱姆法则

令$M$ = $\begin{bmatrix}$
-D & (B-A) & (C-A)
\end{bmatrix}, 也就是将上式的$O-A$作为替换项，依次替换$M$中的$-D$，$B-A$, $C-A$，得到$M_t, M_u, M_v$。那么最终$t = \frac{M_t}{M}, u = \frac{M_u}{M}, v = \frac{M_v}{M}$。理所当然的可以写为下面这类形式

$\left[ \begin{array}{r} t \ u \ v\end{array}\right] = {1 \over \left[ \left| \begin{array}{r} -D & (B-A) & (C-A) \end{array}\right| \right]}$
\left[ \begin{array}{c}
\left| \begin{array}{c} O-A & B-A & C-A\end{array}\right| \
\left| \begin{array}{c} -D & O-A & C-A\end{array}\right| \
\left| \begin{array}{c} -D & B-A & O-A\end{array}\right| \
\end{array}\right]

原文在这里使用了变量替换，本质一致

Scalar triple product很好理解，本质上就是计算一个平行六面体的体积，无论是先进行何方向上的叉积点积，最终计算的都是平行六面体的体积。

纯量三重积最重要的一点就是$A \ldotp (B \times C) = det \left[ \begin{array}{c}$
\begin{array}{c} a_1 & a_2 & a_3\end{array} \
\begin{array}{c} b_1 & b_2 & b_3\end{array} \
\begin{array}{c} c_1 & c_2 & c_3\end{array} \\
\end{array}\right]

$D = B \times C = \left[\begin{array}{c} b_yc_z-b_zc_y & b_zc_x-b_xc_z & b_xc_y-b_yc_x\end{array}\right]$

那么

$A \ldotp D = \left[ \begin{array}{c} a_x& a_y & a_z \end{array} \right] \ldotp D = a_xb_yc_z + a_yb_zc_x + a_zb_xc_y - a_zb_yc_x - a_yb_xc_z - a_xb_zc_y$

对比行列式结果与上式，发现一致。

根据纯量三重积的结论，可以推出

$\left[ \begin{array}{r} t \ u \ v\end{array}\right] =$
{1 \over {(D \times (C-A)) \cdot (B-A)}}
\left[ \begin{array}{c}
((O-A) \times (B-A)) \cdot (C-A) \\
(D \times (C-A)) \cdot (O-A) \\
((O-A) \times (B-A)) \cdot D
\end{array}\right]

注意这里已经通过交换叉积顺序抵消$-D$前方的负号了

---

最终实现代码

bool rayTriangleIntersect( 
    const Vec3f &orig, const Vec3f &dir, 
    const Vec3f &v0, const Vec3f &v1, const Vec3f &v2, 
    float &t, float &u, float &v) 
{ 
    // ---1
    Vec3f v0v1 = v1 - v0; 
    Vec3f v0v2 = v2 - v0; 
    Vec3f pvec = dir.crossProduct(v0v2); 
    float det = v0v1.dotProduct(pvec); 

    // if the determinant is negative the triangle is backfacing
    // if the determinant is close to 0, the ray misses the triangle
    if (det < kEpsilon) return false; 
    float invDet = 1 / det; 

    // ---2
    Vec3f tvec = orig - v0; 
    u = tvec.dotProduct(pvec) * invDet; 
    if (u < 0 || u > 1) return false; 

    // ---3
    Vec3f qvec = tvec.crossProduct(v0v1); 
    v = dir.dotProduct(qvec) * invDet; 
    if (v < 0 || u + v > 1) return false; 

    // ---4
    t = v0v2.dotProduct(qvec) * invDet; 

    return true; 
#else 
    ... 
#endif 
} 

这里第一步先行计算${1 \over {(D \times (C-A)) \cdot (B-A)}}$中的分母是否为0，来判断三角形平面与射线方向是否平行。从分母的几何意义上说，根据交换律本质上原式也等价于$D \ldotp ((C-A)\times(B-A))$，而D的右项即为法线(具体看左右手系)

我这里选取右手系，因为计算得到的法线为反的，因此这里在判断$det$的正负时要看具体左右手系。

得到了det之后就按照克莱姆法则依次计算出t, u, v即可。先行判断$u, v$的取值范围就是用于判断相交点P是否在三角形内。